function [h_dual_metric, riemann_metric] = riemann_metric(Y, laplacian)
    n_samples = size(laplacian, 1);
    n_dim_Y = size(Y, 2);
    h_dual_metric_full = zeros(n_samples, n_dim_Y, n_dim_Y);
    n_use = n_dim_Y;
    
    for i = 1: n_dim_Y
        for j = i: n_dim_Y
            yij = Y(:, i) .* Y(:, j);
            h_dual_metric_full(:, i, j) = 0.5 * (laplacian * yij - Y(:, j) .* (laplacian * Y(:, i)) - Y(:, i) .* (laplacian * Y(:, j)));
        end
    end
    
    for j  = 1: n_dim_Y - 1
        for i = j + 1: n_dim_Y
            h_dual_metric_full(:, i, j) = h_dual_metric_full(:, j, i);
        end
    end

    riemann_metric = h_dual_metric_full;
    for i = 1: n_samples
        riemann_metric(i, :, :) = pinv(squeeze(h_dual_metric_full(i, :, :)));
    end
    
    riemann_metric = riemann_metric(:, 1: n_use, 1: n_use);
    h_dual_metric = h_dual_metric_full(:, 1: n_use, 1: n_use);
end
    
